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Gang  Bao2  and  William  W.  Symes 
Department  of  Mathematical  Sciences 
Rice  University 
Houston,  Texas  77251-1892 


Abstract 

A  simple  algorithm  is  described  for  computing  general  pseudo-differential  operator 
actions.  Our  approach  is  based  on  the  asymptotic  expansion  of  the  symbol  together 
with  the  Fast  Fourier  Transform  (FFT).  The  idea  is  motivated  by  the  characterization 
of  pseudo- differential  operator  algebra.  We  show  that  the  algorithm  is  efficient  through 
analyzing  its  complexity.  Some  of  numerical  experiments  are  also  presented. 

Key  words,  pseudo-differential  operators,  fast  Fourier  transform,  spatially  varying 
filters,  microlocal  cut-off,  data  processing 

Abbreviated  title:  COMPUTATION  OF  PSEUDODIFFERENTIAL 
OPERATORS 

AMS(MOS)  subject  classification:  35S05,  65T20,  86A22 


1  Introduction 

The  theory  of  pseudo-differential  operators  (ty.D.O.s)  has  made  many  important  contri¬ 
butions  to  the  development  of  partial  differential  equations.  It  provides  a  natural  way  to 

1This  work  was  partially  supported  by  the  National  Science  Foundation  under  grant  DMS  89-05878,  by 
the  Office  of  Naval  Research  under  contracts  N00014-J-89-1 115,  by  AFOSR  89-0363,  and  by  the  Geophysical 
Parallel  Computation  Project  (State  of  Texas). 

2current  address:  Institute  for  Mathematics  and  Its  Applications,  University  of  Minnesota,  514  Vincent 
Hall,  Minneapolis,  MN  55455 
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decompose  a  differential  operator  which  may  be  difficult  to  study  directly  into  several  pieces 
with  simple  structure.  A  precise  way  to  describe  propagation  of  singularities  for  differential 
equations  is  in  terms  of  ty.D.O.s.  Pseudo-differential  operators  may  be  viewed  as  spatially 
varying  filters  with  simple  asymptotics  at  high  frequencies.  Pseudo-differential  operators 
differentiate  waves  and  wave-like  signals  according  to  directions  of  propagation.  Pseudo¬ 
differential  operators  also  arise  naturally  in  diverse  fields  (often  under  different  names!)  such 
as  Wave  Propagation,  Electrical  Engineering,  and  Geophysics. 

Although  the  theory  of  pseudo-differential  operators,  or  more  generally  microlocal  anal¬ 
ysis,  has  been  well  established  since  60’s,  little  attention  appears  to  have  been  paid  to  the 
computation  of  pseudo-differential  operators.  In  this  paper,  we  present  a  simple  algorithm 
for  computation  of  general  D.O .  actions.  Our  idea  is  based  on  the  following  characteriza¬ 
tion  of  'FD.O.s: 

Fact:  The  pseudo-differential  operator  algebra  is  generated  by  all  differential  operators  and 
all  powers  of  the  Laplacian. 

More  precisely,  partial  differential  operators  and  many  functions  of  these  (inverse,  powers, 
•  •  •)  are  included  in  ty.D.O.s  in  the  high  frequency  asymptotic  sense. 

See  Kohn  &  Nirenberg  [3]  for  a  detailed  discussion. 


2  Pseudo- differential  Operators 

Here,  we  shall  give  a  brief  introduction  to  a  class  of  'll. D.O. s.  For  a  complete  account 
of  ty.D.O.s  as  well  as  the  calculus,  the  reader  is  referred  to  Taylor  [5],  Nirenberg  [4],  or 
Hormander  [2]. 

We  begin  with  the  introduction  of  Fourier  transform  and  inverse  Fourier  transform.  The 
Fourier  transform  acting  on  a  “nice”  function  u  defined  in  1R"  is: 

J-(u )  =  u  =  (27r)-n  J  u(x)e~iX'^dx 

and  the  inverse  Fourier  transform  is  defined  by 

r'w  =  J  . 
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Remark.  A  number  of  algorithms  for  numerical  computation  of  (discrete)  Fourier  transforms 
have  been  made  available.  We  shall  use  a  version  of  the  Fast  Fourier  Transform  in  our 
numerical  work.  A  detailed  description  may  be  found  in  Conte  and  De  Boor  [1].  See  also 
Van  Loan  [6]  for  the  most  recent  development  in  the  field. 

Pseudo-differential  operators  are  usually  defined  in  terms  of  symbols,  which  are  smooth 
functions  of  both  space  and  frequency  variables  satisfying  certain  estimates.  More  precisely, 
q(x,()  is  a  member  of  the  symbol  class  S'™0(IRn)  iff  q(x,  £)  is  a  smooth  function  and  for  any 
compact  subset  K  of  IR”,  and  real  a,  /3,  there  exists  a  constant  CV, <*,/?,  such  that 

|b;d?«(*,oi  <  cx,aA  1  + 


for  all  x  €  K  and  £  €  IRn. 

In  this  paper,  we  shall  confine  ourselves  to  a  subclass  of  £™0,  the  class  Sm ,  which  is 
the  most  natural  class  and  sufficient  for  many  applications.  A  function  q(x,()  is  in  Sm  if 
q(xiO  G  S™o  and  there  are  smooth  qm-j{x,£),  homogeneous  of  degree  m  —  j  in  £  for  |£|  >  1, 
t.e., 

gm-i(*,r£)  =  rm_Jgm-j(z,£),  If  I  >  1  ,  r  >  1 


such  that 


in  the  sense  that 


j>  0 

q(x,0  -  £}qm-j(x,0  e  5,1m0-Ar-1 , 

1=0 


(2.1) 


where  gm(x,£)  is  called  the  principal  symbol,  or  principal  part  of  q(x,£)  which  carries  the 
most  important  information  about  q. 

Then  the  operator  Q  defined  by 


Q(x,  Dx)u  =  J  q(x,Z)u(t)e,x<d( 

is  called  a  ty.D.O.  of  order  morQg  OP5'm(lRn). 

In  particular,  differential  operators  with  smooth  coefficients  are  ty.D.O.s.  Indeed,  for 
such  a  differential  operator  of  order  m,  the  corresponding  symbol  is  a  polynomial  in  £  of 
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degree  m,  and  consequently  is  a  symbol  in  Sm.  The  asymptotic  expansion  of  a  symbol,  (2.1), 
is  unique  up  to  smoothing  operators. 

3  Algorithm 

In  this  section  we  describe  the  algorithm  explicitly.  Its  complexity  will  be  examined  in  the 
section  that  follows.  For  the  sake  of  simplicity,  we  shall  only  describe  the  idea  of  computing 
two-dimensional  ty.D.O. s.  Some  obvious  modifications  may  be  made  to  compute  ty.D.O.s 
of  arbitrary  dimension.  Throughout,  we  shall  always  assume  that  the  action  of  Q  on  u  is 
meaningful.  The  precise  conditions  may  be  found  in  any  one  of  the  above  references. 

Given  the  symbol  q(x,  z,  £,  q)  of  a  ^.D.O.  Q(x,  z,  Dx ,  Dz)  6  OPSm  and  a  function  u(x,z). 
Let  us  assume  that  the  asymptotic  expansion  of  the  symbol  is  given  by 

<l(x,  z,  n)  ~  Yj  2,  L  V)  (3-2) 

j>o 

where  qm-j{x,z,r(,rq)  =  rm~3 qm_j{x,z,£, q)  for  (|£|,  |i/|  >  1).  Again,  qm  denotes  the 
principal  symbol  of  q. 

Knowing  the  asymptotic  expansion  of  q  we  compute  the  ty.D.O.  action  Qu  through 
computing  Qm-jU  for  j  >  0.  We  describe  the  calculation  for  j  =  0,  as  it  is  representative. 
That  is  we  will  describe  an  algorithm  for  computing  the  action  of  the  principal  part  of  a 
ty.D.O..  Evidently  this  algorithm  could  be  applied  recursively  to  compute  general  ty.D.O. 
actions.  However  the  principal  part  gives  the  dominant  effect  on  high  frequency  inputs, 
which  is  most  important  for  our  intended  applications.  Thus  for  us,  computation  of  the 
principal  part  above  is  sufficient. 

Let  ^  =  wcosO,  T]  =  usind,  lo  =  \/£2  +  rj2 ,  then  the  homogeneity  of  qm  in  £  and  q  yields 
that 

qm(x,z,t,i ])  =  qm(x,z,ucos9,usin6 )  =  umqm(x,z,0 ) 
where  qm(x,z,6)  is  defined  to  be  qm(x,z,cos0,sin6). 

Since  qm  is  periodic  in  0 ,  the  Fourier  series  expansion  can  be  employed  to  give 

oo  K  j  2 

qm(x,z,0)=  ci(x,z)e,ie  ~  £  q(x,  z)etW 

l=- oo  l=-K/ 2 
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(3.3) 


K/2 

=  ct(x,z)(cosQ  +  isinO)1  . 

l=-K/2 

It  follows  that  from  the  definition 


Qmu  ~  f  I  d(drje'^+z^  ^  u>mcl(x,z)(cosO  +  isin9)lu(^q) 
J  J  l=-K/2 

Kl 2  ,  , 

=  ^2  ci(x,z )  /  I  d£dr]umet(xt+ZT,\cos6 isin9)lu(£,T)) 

t  _  ry  trx  J  J 


l=-K/2 
K/2 

l=-K/2 


where  to  obtain  the  last  equality,  we  have  used  the  relations  cosO  =  (/to  and  sin9  —  tj/u ;. 
Observe  that  um~l  is  the  symbol  of  the  (m  —  /)/2- power  of  the  (negative)  Laplacian,  while 
£  and  Tj  are  symbols  of  differential  operators  Dx  =  —idx  and  Dz  —  —idz,  respectively. 

The  procedure  implicit  in  the  above  formulae  leads  to  an  algorithm  to  evaluate  Qmu 
approximately,  as  follows.  Assume  that  u  is  sampled  on  a  discrete  grid, 


Uij  =  u(x0  +  (i  -  l)Ax,  zo  +  (j  -  1) A z) 
i  =  1,-  •  -,M,  j  =  1,-  •  *,JV 


with  spacings  Ax,  A z  >  0.  Assume  similarly  that  a  sampling  of  qm  is  given, 


Qi,j,k  =  qm{x o  +  (t  -  1)  Ax,  z0  +  (i  -  l)Az,  kA9) 
i  =  1,***,M,  j  =  1,  •  •  -,N,  k  =  -K/2,---,I</2 

with  X  =  ( M  —  l)Ax,  Z  —  (N  —  l)Az,  the  sample  rates  in  the  frequency  domain  are 
A£  =  1/X,  At)  =  1/Z,  so  the  (unaliased)  samples  of  the  symbols  of  the  Laplacian,  DX)  and 
Dz  are 

Op,r  =  2'K^'pAO,'1  +  (rArj)2 
=  27rpA£ 

ZPtr  =  2TcrArj 

p  =  -M/2,  •  •  -,  M/2,  r  =  —N/2,  ■  •  •,  N/2 
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respectively. 


Procedure  for  Computing  Qmu 

1.  Compute  the  discrete  Fourier  transform  U  of  U. 

2.  For  each  i  €  { 1 ,  -  -  *,M}  and  j  6  {1,*  •  • ,N },  compute  the  discrete  Fourier  transform 

QiJ  =  { Q i,j,l } l=—K/2  Qij  =  {Qi,i,k}k=-K/2- 

3.  Initialize  (QU)ij  =  0.0,  i  —  1,  •  •  •,  M,  j  =  1,  •  •  •,  N. 

DO  /  =  -K/ 2,  K/ 2 

a  compute  the  inverse  transform  of 

n;;‘(EPir  +  iz„?u„ 

for  p  =  -M/2,  •  •  -,  M/2  and  r  =  -N/ 2,  •  •  •,  AT/2. 
b  accumulate 

(QU)ij  =  (QU)ij  + 

END  DO 


4  Complexity  Analysis 

We  return  to  the  general  case.  The  complexity  will  be  analyzed  by  the  number  of  multipli¬ 
cations.  We  also  make  a  few  remarks  about  accuracy  of  the  algorithm. 

The  direct  method  of  computing  the  ty.D.O.  action  is  by  straightforward  discretization 
of  the  definition. 


Qmu  =  J  j  qm(x,0u(0eix<dZ 


is 


Let  us  assume  that  the  input  function  is  discretized  on  a  regular  d-dimensional  grid,  as 
the  symbol  qm.  We  denote  by  N  the  number  of  grid  points  in  each  direction,  assuming 
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these  are  roughly  similar.  Assuming  also  that  the  discrete  Fourier  transforms  are  computed 
using  a  Fast  Fourier  transform  algorithm,  we  then  have  the  following  result. 

Lemma  4.1  The  direct  algorithm  has  0(N2dlogN )  complexity. 

This  is  an  immediate  consequence  of  the  well  known  fact  that  the  FFT  exhibits  O(NlogN) 
complexity,  where  N  is  the  length  of  the  input  sequence. 

We  next  discuss  the  complexity  of  the  new  algorithm.  For  simplicity,  we  once  again 
consider  the  2 -D  case.  The  approximate  complexity  orders  of  the  steps  in  the  algorithm 
proposed  above  are: 

1.  N2logN 

2.  N2KlogK 

3.  a.  KN2logN ;  b.  KN2. 

Hence  the  total  complexity  is  0(KN2(logN  +  logK ))  in  two  dimension. 

In  general,  when  the  number  of  dimension  is  d,  a  similar  calculation  will  yield 

Lemma  4.2  The  new  algorithm  exhibits  Nd(logN  +  logK))  complexity. 

Remarks.  The  new  algorithm  is  significantly  superior  to  the  direct  method.  The  number 
of  terms  K  in  the  finite  0-Fourier  series  approximation  of  qm  ought  to  be  chosen  so  that 
the  sup  norm  error  is  small.  Then  the  error  in  the  computation  of  Qm,  modulo  compact 
operators,  will  also  be  small  (Taylor  [5],  p.  52).  In  particular  the  error  will  be  small  for 
oscillatory  inputs  u.  Note  that  K  is  completely  independent  of  N  in  this  regard.  Thus  in 
effect  the  complexity  of  our  algorithm  is  0(NdlogN)\ 

5  Numerical  Experiments 

In  this  section,  we  present  the  results  of  some  numerical  experiments  carried  out  with  the 
ty.D.O.  calculation  algorithm.  We  calculated  actions  involving  several  microlocal  cutoff 
operators,  i.e.  operators  whose  symbols  are  asymptotically  1  in  some  conic  set  and  asymp¬ 
totically  zero  in  the  complement  of  another  conic  set  (essential  support  or  aperture).  These 
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are  the  simplest  indecomposable  order  zero  \Er. D.O.s .  The  examples  exhibit  some  of  the 
interesting  features  of  Vt. D.O.s. 

We  began  with  convolutional  ty. D.O.s,  which  are  ty. D.O.s  that  are  independent  of  spatial 
variables.  Obviously  convolutional  operators  are  natural  extension  of  differential  operators 
with  constant  coefficients.  For  this  class  of  operators,  it  is  easy  to  show  that 

HQu)  =  Q((HO , 

which  is  useful  in  verifying  the  code.  In  fact,  according  to  this  simple  identity,  one  can  recover 
the  symbol  from  Qu  and  it.  A  symbol  that  characterizes  a  microlocal  cutoff  is  specified  by 
Figure  1,  where  the  symbol  was  designed  to  be  a  C2  function.  The  next  figure  displays 
the  symbol  function  in  terms  of  the  angle  0.  From  this  one  dimensional  array,  the  fy.D.O. 
algorithm  was  employed  to  compute  the  action,  and  hence  the  2 -D  symbol  function  Q.  The 
result,  Figure  3,  shows  the  symbol  when  the  number  of  terms  in  Fourier  series  expansion 
of  the  symbol  k  =  4.  It  is  easy  to  see  that  the  symbol  in  Figure  3  illustrates  the  right 
direction  but  wrong  amplitude  within  the  aperture.  As  we  increased  k,  the  recovery  of  the 
symbol  became  better  and  better.  Figure  4  shows  that  the  symbol  was  perfectly  recovered 
after  several  steps.  Again,  we  want  to  emphasize  that  the  number  k  only  depends  on  the 
smoothness  of  the  symbol. 

The  next  experiment  concerned  the  rotation  of  apertures  for  convolutional  operators. 
The  function  plotted  in  Figure  5  is  a  slightly  smoothed  characteristic  function  of  a  circle. 
We  applied  a  ty.D.O.  cutoff  whose  symbol  was  given  in  Figure  6  to  this  function.  Just 
as  the  theory  predicts,  the  high  frequency  information  of  the  resulting  function  (Figure 
7)  was  preserved  within  the  aperture.  We  then  rotated  the  symbol  (Figure  8  and  Figure 
10),  and  again  the  high  frequency  information  were  preserved  in  Figure  9  and  Figure  11, 
respectively.  These  examples  are  only  illustrative,  as  the  discrete  Fourier  transform  allows 
a  very  simple  and  fast  computation  of  convolution  operators.  Our  next  example  is  meant 
to  illustrate  the  success  of  our  algorithm  with  nonconvolutional  .D.O.s.  Figure  12  shows 
the  symbol  of  a  2-D  ty.D.O.  which  is  spatially  varying  (in  z  direction).  It  was  generated 
from  q(z,6)  =  q0(6  +  86sin(nz/ zmax))  where  q0  was  given  in  Figure  6,  86  was  selected  to  be 
7r/2,  and  z  €  [0,  zmax].  Thus,  as  z  increases,  the  symbol  rotates  smoothly,  in  particular  the 
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symbol  will  be  equal  to  q0  as  2  reaches  its  maximum.  Once  again,  we  used  the  function  u 
in  Figure  5.  The  result,  as  shown  in  Figure  13,  agreed  with  the  theory.  Observe  that  the 
aperture  is  vertical  for  z  near  0.  Then  the  symbol  rotates  as  z  increases,  so  we  start  to  see 
some  high  frequency  horizontal  components.  When  z  is  getting  close  to  its  maximum,  the 
symbol  rotates  back,  and  the  aperture  becomes  vertical  again. 

Our  final  example  demonstrates  an  important  application  of  the  ty.D.O.  algorithm  to  the 
seismic  data  processing  in  reflection  seismology.  The  basic  objective  of  all  seismic  processing 
is  to  convert  the  information  recorded  in  a  field  into  a  form  that  most  greatly  facilitates 
geological  interpretation  of  such  field.  Evidently,  real  reflection  data,  which  carry  most 
of  the  information  of  the  mechanical  properties  of  the  earth,  are  what  geophysicists  most 
interested  in  obtaining  through  this  process.  Thus,  an  essential  object  of  the  processing 
is  to  eliminate  or  suppress  all  signals  not  associated  with  reflections.  Figure  14  displays  a 
seismogram,  i.e.  the  recorded  seismic  data  at  receivers  on  the  surface  of  the  earth  after  an 
energy  source  is  fired.  As  one  can  see  clearly  that  the  dark  area  that  contains  very  strong 
signals  represents  the  early  arrivals  (direct  and  head  waves).  In  the  area  below,  there  are 
other  signals  (reflections)  which  are  not  nearly  as  strong  as  the  early  arrivals.  Unfortunately, 
the  direct  and  head  waves  do  not  penetrate  the  earth,  so  contain  no  information  about 
the  subsurface  that  we  are  interested  in.  The  useful  reflection  energy  is  only  contained  in 
the  lower  part.  This  can  be  observed  more  clearly  if  one  increase  the  amplitude  of  the 
seismogram  as  in  Figure  15.  The  question  arose:  Can  one  remove  the  early  arrivals  and  yet 
keep  the  useful  information  of  reflections?  Applying  the  tf.D.O.  computation  algorithm,  we 
were  able  to  design  a  microlocal  cutoff  (ty.D.O.)  whose  action  on  the  seismogram  is  shown 
in  Figure  16.  The  result  appears  to  be  very  encouraging.  The  amplitude  of  the  early  arrivals 
is  reduced  dramatically,  and  meanwhile  information  of  reflections  is  well  preserved.  We  tried 
the  same  ty.D.O.  filter  on  the  resulted  data  set  one  more  time  to  obtain  Figure  17.  The  early 
arrivals  are  essentially  gone,  while  again  most  of  the  reflections  are  preserved.  We  believe 
the  noise  left  in  the  region  where  the  early  arrivals  resided  was  caused  by  numerical  scales, 
hence  can  be  eliminated.  This  processing  technique  is  actually  used  in  reflection  seismology, 
where  it  is  called  uf-k  dip  filtering”,  e.g.  Yilmaz  [7],  pp  69-78.  Our  ty.D.O.  algorithm  yields 
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an  accurate  and  efficient  means  of  “spatially  variable  dip  filtering”,  for  which  we  envision 


numerous  uses. 


6  Concluding  Remarks 

A  simple  algorithm  for  the  computation  of  a  class  of  ^.D.O.s  is  introduced  in  this  work. 
We  exhibit  some  of  the  features  of  the  algorithm.  The  complexity  analysis  indicates  that 
the  algorithm  is  much  more  efficient  than  the  direct  computation.  Because  of  the  simple 
structure,  various  massive  parallel  computers  may  be  used  to  implement  this  algorithm 
so  long  as  a  fast  FFT  routine  and  fast  array  operations  are  available.  In  fact,  some  of 
our  numerical  experiments  reported  in  this  paper  were  obtained  by  using  the  Connection 
Machine. 

We  anticipate  many  applications  of  this  algorithm.  For  example,  ty.D.O.s  are  expected 
to  play  an  important  role  in  regularizing  a  class  of  ill-posed  problems  in  multi-dimensional 
wave  propagation  arisen  naturally  in  seismic  inversion,  oil  and  gas  exploration,  and  many 
other  related  geophysical  problems.  Our  experiment  indicates  the  usefulness  of  microlocal 
(or  ty.D.O.)  cutoff  in  seismic  data  processing,  i.e.,  sorting  of  waves  according  to  direction  in 
seismic  data. 

Mathematically,  this  algorithm  should  provide  a  way  to  compute  the  so  called  microlocal 
norms  of  microlocal  Sobolev  spaces,  which  in  turn  would  help  us  testing  the  sharpness 
of  various  results  on  propagation  of  singularities  for  partial  differential  equations.  This 
algorithm  should  also  have  some  impact  on  signal  processing,  where  ty.D.O.s  form  a  class  of 
spatially  varying  filters. 
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Figure  6:  symbol  q0 ■ 


120 


20  40  60  MO  100 

Figure  8:  symbol  rotation  1. 


Figure  10:  symbol  rotation  2. 
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Figure  12:  a  spatially  varying  symbol. 
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Figure  13:  the  action  for  the  spatially  varying  symbol. 


Figure  14:  seismogram. 


Figure  15:  seismogram  with  increased  amplitude. 
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Figure  16:  the  result  after  applying  the  filter. 


Figure  17:  the  result  after  applying  the  filter  twice. 


